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ABSTRACT 

A significant fraction of massive stars are moving supersonically through the interstellar medium 
(ISM), either due to disruption of a binary system or ejection from their parent star cluster. The 
interaction of their wind with the ISM produces a bow shock. In late evolutionary stages these stars 
may undergo rapid transitions from red to blue and vice versa on the Hertzsprung-Russell diagram, 
with accompanying rapid changes to their stellar winds and bow shocks. Recent 3D simulations of the 
bow shock produced by the nearby runaway red supergiant (RSG) Betelgeuse, under the assumption 
of a constant wind, indicate that the bow shock is very young (< 30 000 years old), hence Betelgeuse 
may have only recently become a RSG. To test this possibility, we have calculated stellar evolution 
models for single stars which match the observed properties of Betelgeuse in the RSG phase. The 
resulting evolving stellar wind is incorporated into 2D hydrodynamic simulations in which we model 
a runaway blue supergiant (BSG) as it undergoes the transition to a RSG near the end of its life. 
We find that the collapsing BSG wind bubble induces a bow shock-shaped inner shell around the 
RSG wind that resembles Betelgeuse's bow shock, and has a similar mass. Surrounding this is the 
larger-scale retreating bow shock generated by the now defunct BSG wind's interaction with the ISM. 
We suggest that this outer shell could explain the bar feature located (at least in projection) just in 
front of Betelgeuse's bow shock. 

Subject headings: circumstellar matter — hydrodynamics — stars: evolution — stars: winds, outflows 
— stars: individual: Betelgeuse 



1. INTRODUCTION 

Winds from massive stars have a strong effect on 
their circumstellar environmen t by driving shocks and 
producing parsec- scale shells (|Dvson fc de Vriesl Il972t 
iCastor et al.lll975l ). When the star is moving these shells 
become asymmetric and, if the star is moving super- 
sonically with respect to the interstellar medium (ISM), 
a bow s hock is formed whose shape was calculated by 
iWilkinl (|1996f ) . A significant fraction of massive O-stars 
become runaway stars because of dyn amical ejection 
from clusters and binary disruption (e.g. iGies fc Boltonj 
19861; iPflamm-Altenburg fc Kroupall2006t lEldridge et al.l 
20111 ). Bow shocks around massive stars have been 



found in optical images (jGull fc Sofial Il979h and with 
the IRAS satellite in the mid- and far-infrared (IR) 
(|van Buren fc McCravl[T98l Ivan Buren et al.lll995l ). IR 
observations are very effi cient for detecting and charac- 
terising bow shocks (e.g. iGvaram adzc & Bomansl 120081 : 
iPovich et a ll 120081: iWinston et al.ll2012D. ~~~ 

IRAS observations (jNoriega-Crespo et al.lll997f l of the 
red supergiant (RSG) Betelgeuse (a Orionis) detected 
a bow shock and a linear "bar" feature upstream from 
the bow shock and perpendicular to the star's proper 
motion, show n here in Fig. [J Higher resolution AKARI 
observations (jUeta et al.ll2008f) confirmed these features . 
Recent Herschel observations (|Cox et al.l I2012T) showed 
layered substructure in the bow shock. iLe Bertre et all 
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Figure 1. Betelgeu se's bow shock a nd "bar" at 60/im, from the 
IRAS Galaxy Atlas l lCao et al.H1997Tl . 



<|2012f ) found a thin shell of far- ultraviolet emission at the 
same position as the IR bow shock. They also presented 
VLA 21cm observations with evidence of some atomic 
hydrogen near the bow shock, and stronger emission from 
an inner shell three times closer to Betelgeuse. 
Thin shells and bow shocks are often subject to insta- 
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bilities for realistic stellar winds and ISM parameters, 
requiring m ulti-dimensional simulation s for quantitative 
studies (e.g. iGarcia-Segura et al.|["l996D . Runaway stars 
with constant winds have been ex tensively modelled for 
both fast winds from hot stars dBrighenti fc D'Ercold 
11995b! : IRaga et all H997t iComeron fe Kaoerl Il998^~and 
slow winds from cool stars (e.g. IVillaver et al.1 l2012t 
ICox et all 120121 and references therein) . Bow shocks 
produced by stars evolv ing from red to blue have 
been studied previou sly (jBrighenti fc D'Erc olc 1995a; 
Ivan Marie et al.ll2006l ). but there have been no detailed 
st udies of the opp os ite tr ansition from blue to red. 

iMohamed et al.l (|2012l ) modelled Betelgeuse's bow 
shock using 3D SPH simulations assuming a constant 
stellar wind. By comparing the observed mass (M s ~ 
0.0033 M from AKARI data) and shape of the bow 
shock with simulations they found that the bow shock 
must be very young (< 30 kyr) and not yet in a steady 
state. This implies that either Betelgeuse has recently 
moved into a region of the ISM with properties very dif- 
ferent to where it came from, or that the wind itself has 
changed significantly. If Betelgeuse had recently under- 
gone a transition from a blue supergiant (BSG) or main 
sequence (MS) star to a RSG, then the bar feature could 
be a r emnant from the MS or BSG wind (|Mohamed et al.1 
I2012D . In this paper we test this idea and investigate such 
a transition with 2D hydrodynamical simulations of an 
evolving wind. 

2. DESCRIPTION OF CODE AND METHODS 
2.1. Stellar evolution models 

We consider a 15 M© solar metallicity star whose 
evolution on the Hertzsprung-Russell diagram is shown 
in Fig. EfA). The tip of the RS G branch has prop- 
erties consistent w ith Betelgeuse (jNeilson et al J 12011 at 
ISmith et al.l 120091 and r eferences therein ) . Th e model 
is computed using the lYoon fc Langerl (|2005|) stellar 
evolution code including mass loss but no convective 
core overshooti n g. M ass loss is computed using the 
iKudritzki et al.l (119 89) approximation for hot stars and 
the lde Jager et al.l ()1988D relation at cooler effective tem- 
peratures. Wind vel ocities are approximated using the 
lEldridge et al] (|2006fl prescription, v 2 = /3 w (T)v% sc . The 
value of of f3 w for T < 3600 K was reduced from their 
value of j3 w — 0.125 to f3 w — 0.04 to match the observed 
wind velocity of Betelgeuse (v w ~ 17 kms -1 ). The wind 
velocity, mass- loss rate M, and wind density (M/v w ), 
are plotted in Fig. H^B) for the last ~ 75 kyr of evolu- 
tion. This corresponds to the end of the blue loop and 
evolution along the RSG branch, highlighted in blue in 
Fig. HA). 

2.2. Simulation setup 

Th e hydrodynamics code we have used ()Mackev fc Liml 
2010) is a finite-volume, uniform-grid code; here the 
equations of inviscid, compressible hydrodynamics are 
solved on an axisymmetric 2D grid in (z, R), using a 
second-order-accurate integration scheme (jFalld fl991). 
Microphysical cooling processes provide a source term 
to the energy equation and are included in the code by 
operator splitting. Cooling rates are calculated using the 
collisional ionization e quilibrium (CIE) cooling tables of 
iWiersma et al.1 (|2009f) for solar abundances. For numer- 
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Figure 2. (A) Hertzsprung-Russell diagram for the 15 Mq stellar 
model. The star symbol identifies t = 11.4 Myr, midway through 
the extended blue loop and the starting point of the simulations in 
this paper. The blue part of the curve highlights the final 75 kyr of 
evolution, shown in more detail in panel (B). Open circles identify 
the times of the four snapshots in Fig. \3\ (B) The wind mass-loss 
rate, M in 10 -7 Mq yr -1 , wind velocity, v w in kms -1 , and M/v m 
(proportional to wind density) are shown. 

ical stability the gas temperature is limited to T > 30 
K. 

The stellar wind bound ary c ondition is implemented 
as in iFrever et all (|2003l ) and Ivan Marie et all (|2006f l 
by imposing a freely expanding wind within a 20 grid- 
zone radius of the origin. Wind parameters are up- 
dated from the stellar evolution model every timestep. 
The initial ISM density and pressure were set to po — 
4.676 x 10 -25 gcm~ 3 and p = 4.676 x 10~ 13 dyne cm -2 , 
respectively, corresponding to a H number-density of 
«h = 0.2 cm -3 for hh = p/2.338 x 10~ 24 g and a temper- 
ature T — 7700 K for ionised gas containing 10% helium 
by number. The star's velocity through the ISM was set 
to = 50 kms -1 . 

The BSG phase has an increasing mass-loss rate of 
M ~ (2 — ► 3) x 10~ 7 M Q yr~ 1 as it approaches the 
transition to RSG, and a decreasing terminal velocity 
of v w ~ 400 -> 300 kms" 1 (see Fig. EJ. Equating the 
wind and ISM ra m pressures gi ves the well-known stand- 
off distance (e.g. I Wilkin! [19961) : i? S o - 0.6 - 0.7 pc for 
the BSG wind just quoted. The simulation domain was 
set to z e [-2.56,1.28] pc and R 6 [0,1.92] pc, using a 
coordinate system where the star is at the origin and the 
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upstream direction is +z; the code works in the reference 
frame where the star is stationary. The upstream bound- 
ary condition was inflow, the symmetry axis was reflec- 
tive, and the other boundary conditions were outflow- 
only. The grid contained 1536 x 768 zones, giving a zone 
size Sx — 0.0025 pc and a wind-boundary radius of 0.05 
pc. A lower resolution simulation verified the numerical 
convergence of our results. The simulation was started at 
t = 11.4 Myr in the evolution model, about 400 kyr be- 
fore the transition from BSG to RSG begins (see Fig. [5]). 
The ISM advects across the domain in 75 kyr, so 400 
kyr is more than sufficient for any effects of the initial 
conditions to be erased. 

3. RESULTS 

Snapshots from the simulation showing nu (above) and 
T (below) are presented in Fig. [3] at the evolutionary 
times indicated. A more quantitative plot of the evolving 
density field along the symmetry axis is shown in Fig. 01 
where the four curves labelled A-D correspond to panels 
A-D of Fig. H The first snaphot (A) shows the CSM 
at the end of the BSG evolutionary phase, with a sin- 
gle bow shock comprising the forward shock at r = 0.95 
pc, cooling-regulated compression in the shocked ISM at 
r ~ 0.9 pc, the contact discontinuity separating wind 
material from the ISM at r = 0.82 pc, and the re- 
verse/termination shock in the BSG wind at r = 0.63 
pc. The forward shock is stable, with a density jump 
of ~ 4x across the shock, and a subsequent increase in 
density as post-shock gas cools. Kelvin-Helmholtz insta- 
bilities are excited by shear across the contact disconti- 
nuity. Further from the symmetry axis these instabilities 
can significantly affect the bow shock position and cur- 
vature. 

During the BSG phase the shocked wind is approxi- 
mately adiabatic, but the shocked ISM lies between the 
adiabatic and isothermal limits for the cooling function 
and ISM density considered. The strong contact discon- 
tinuity ensures that the shell thickness is set only by the 
cooling of the shocked ISM layer. Without cooling the 
shell is w 3x thicker; with enhanced cooling it is thin- 
ner and ve ry unstable. Our mod el appears to use similar 
cooling to Wa reing et al.l (120071) b ut possibly somewhat 
weaker than IVillaver et al.l (|2012T ). The uncertain bal- 
ance between heating and cooling at T < 10 4 K at these 
low densities makes the shell thickness difficult to model 
reliably; additionally magnetic pressure could also play 
a role. 

At time (B), the wind density has increased at least 
tenfold out to r ~ 0.15 pc and there is a shell with en- 
hanced density in the wind at r ~ 0.3 pc from the earlier 
local maximum in the wind density at t — 11.868 Myr 
(Fig. 03). The BSG bow shock is no longer supported 
by the fast wind's ram pressure, so it retreats towards 
the slowly expanding RSG wind, preceded by the reverse 
shock. The RSG wind has no shell because v w decreases 
monotonically with time during the transition (so suc- 
cessive layers of wind material don't interact). 

Between times (B) and (C) the hot bubble still has a 
large sound speed, so the reverse shock collapses onto 
the RSG wind rapidly, compressing the shell at r ~ 
0.3 pc, and creating a thinner, denser shell at the edge 
of the RSG wind by 11.8829 Myr. The mass of the in- 
ner shocked-shell grows linearly with time from zero at 



11.882 Myr to 0.002 M Q at 11.892 Myr, when it reaches 
the upstream contact discontinuity and merges with it. 

At time (C) the thin shell surrounding the RSG 
wind is shaped like a bow shock, with mass and phys- 
ical size similar to tha t determined from AKARI data 
(jMohamed et al.l I2012H . A second, much weaker shell 
is also seen at a slightly larger radius. Both the BSG 
forward shock and contact discontinuity are still being 
swept downstream, but the weak shell-like feature at 
z = 0.44 pc in curve C shows that gas pressure in the 
compressed hot bubble is impacting the contact discon- 
tinuity. We identify the retreating contact discontinu- 
ity as a counterpart to the "bar" upstream from Betel- 
geuse's bow shock. This transitory phase, from when the 
RSG wind encounters the BSG reverse shock to when it 
reaches the BSG contact discontinuity, lasts 10 kyr. The 
two shells have merged by t = 11.8939 Myr. 

After time (C) the density structure is that of lay- 
ered shells, which soon merge as the RSG wind pushes 
through the BSG bow shock, and by the end of the star's 
life at time (D) the RSG wind has reached the undis- 
turbed ISM and is forming a new bow shock. The RSG 
wind is also bounded by a shell downstream because of 
the pressure of the collapsing BSG bubble. It would take 
~ 60 kyr (from the time it loses pressure support) for the 
BSG shock to be completely swept downstream on the 
simulation domain shown, so the transition phase is far 
from complete at the pre-SN stage. 

4. DISCUSSION AND COMPARISON TO BETELGEUSE 

IR observations o f Betelgeuse' s CSM 

(iNoriega-Crespo et al.l[l997t lUeta et al.M2008b fCox et al.l 
12011 " show that the standoff distance of the bow shock 
is w 0.3 pc, and that of the bar is * 0.45 pc (for d = 200 
pc). Fig. [3] shows that a BSG bow shock which collapses 
as the star becomes a RSG produces strikingly similar 
structures. In this model the "bow shock" around 
Betelgeuse is actually the collapsing BSG reverse shock 
interacting with the expanding RSG wind, and the bar 
is the BSG contact discontinuity. 

The total gas mass in Betelgeuse's combined wind and 
bow shock is estimate d at 0.033 M from IRAS data 
(jMohamed et al.ll2012t ); in our model at time (C) the 
total mass within r < 0.3 pc is 0.030 M , in good 
agreement. The bo w shock mass is M« ~ 0.0033 M Q 
from AKARI data (jMohamed et alj 120121), or possi- 
bly as small as 3.8 x 1O" 4 M (jCox et al.ll20ll . This 
is > 20 x less massive than steady-state bow shocks 
(M s ~ 0.05 — 0.15 Mf?)) from simulations of runaway 
RSGs (jMohamed et al.l l2012f) . a problem which moti- 
vated the study presented here. In our model the expand- 
ing RSG wind effectively encounters a vacuum initially, 
so there is no swept up shell from the early expansion and 
consequently its mass is much lower than a normal RSG 
bow shock, M s < 0.002 M . Our model therefore pro- 
vides a good explanation for the mass of the bow shock, 
and is consistent with the RSG wind m ass. 

Herschel observations of Betelgeuse (jCox et alJ I2"012h 
show that the bow shock seems to have a layered struc- 
ture. The local maximum in wind density during the 
BSG— 5-RSG transition (Fig. [5]) remains visible in our sim- 
ulation (Figs. 131 and HJ) for some time outside the inner 
shell, and such modulations in the wind properties could 
explain the observed layers. 
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Figure 3. Gas number density njj (upper half-plane, in cm s ) and temperature (lower half-plane, in Kelvin) plotted on logarithmic scales 
at times t = 11.85, 11.875, 11.8875, and 11.9073 Myr from left to right, respectively. Only part of the simulation domain is shown. 
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Figure 4. Gas density (cm -3 ) as a function of position (in pc, 
relative to the star) along the symmetry axis. Curves are labelled 
by the simulation time (in Myr) of the snapshot; each curve is 
offset by +1 dex from the previous one. Labels A-D correspond to 
the snapshots shown in panels A-D of Fig. [3] 
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Figure 5. Projected dust luminosity in the IRAS 65 /im band 
(above) and projected mass density (below) at t = 11.8875 Myr. 
The dust luminosity is shown in units of ergs -1 per pixel, and 
the projected mass in gem -2 , both on a linear scale. 20 contours 
are overlaid using the same linear scale. Not all of the simulation 
domain is shown. 

There are two upstream discontinuities that could cor- 
respond to the bar seen upstream from Betelgeuse's bow 
shock, namely the BSG forward shock and contact dis- 
continuity (or both, if the bow shock were a thin shell). 
For the simulation presented here the contact disconti- 
nuity provides the most likely counterpart, based on its 
position and shape at the time the inner shocked-shell 
is visible. We calculate the observed properties of the 
2-D models (they have cylindrical symmetry) to test this 
hypothesis. The projected dust-emission intensity (up- 
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per half-plane, in ergs 1 per pixel) and projected mass 
(lower half-plane, in gem -2 ) are shown in Fig. [S] at time 
(C) (11.8875 Myr). Dust emission in the IRAS 65 /mi 
waveband is calculated assuming radiative equilibrium, 
and using t he same assumption s for dust size and com- 
position as iNeilson et al.l J2010) with a dust-to-gas ra- 
tio of 1/200. The inner shocked-shell and a prominent 
outer shell are still apparent, both in projected density 
and in dust emission. The r~ 2 decrease in dust emission 
makes the contact discontinuity look a little more like a 
bar (compared to projected density) because the strongly 
curved parts at larger radii are fainter. The total dust 
luminosity of the RSG wind, bow shock, and bar from 
this model are lower than the IRAS/ AKARI observa- 
tions, suggesting radiative excitation from Betelgeuse's 
radiation may not be the dominant process exciting the 
dust emission. Converting from simulated mass density 
to far-infrared emission is complicated, however, because 
the dust composition of the ISM, BSG and RSG winds 
are probably all somewhat different. Increasing the dust- 
to-gas ratio would increase the predicted luminosity. It is 
also not clear what contribution C and O fine-structure 
lines make to the far-IR emission. We have only consid- 
ered Betelgeuse's radiation as the energy source for the 
dust emission, but radiation produced within the bow 
shocks and hot bubble probably also contribute. 

The observed bar has no discernible curvature, but a 
strong prediction of our model is that it must eventually 
begin to curve back into a bow shock shape, if physically 
associated with Betelgeuse. The radius of curvature can 
be larger than the current standoff distance of the bar 
would suggest if it is a relic structure that is currently 
being overtaken by the star (cf. Fig. [3]). Instabilities pro- 
vide another mechanism for reducing bow shock curva- 
ture (compare panels A and B of Fig. [3]): the shock flexes 
periodically as vortices are generated and shed down- 
stream, and the curvature of the contact discontinuity 
can consequently fluctuate quite dra matically. This has 
also been seen in previous work (e.g. [Comeron fc Kaperl 
19981: iWareing et aI1l2007t iVillaver et al.ll2012D . 

The transition timescale, from when the RSG wind 
appears to when the BSG wind's bow shock has dis- 
appeared, can be quite long and is beyond the end of 
the star's life in this model. The BSG bow shock loses 
pressure support on a timescale comparable to the 10-15 
kyr required for the star to expand to a RSG. Subse- 
quently, if the star is moving with w» = 50kms _1 and 
the wind has Rso = 1 pc, the time to cross this distance 
is ~ 20 kyr but it can take three times as long for the bow 
shock to stall and be swept far downstream. For larger 
Rso and lower it is possible for the observable bow 
shock transition to last ~ 100 kyr, a significant fraction 
of the massive star's post-MS evolution, hence it may 
not be unusual to observe a RSG during this transition. 
The transition from MS to RSG (not presented here) is 
qualitatively similar, differing in the wind strengths and 
speed of transition. 

A number of uncertainties affec t our model. Betel- 
geuse's mass is poorly constrained (jNeilson et al.ll20l"Tal) 
and the properties of the proposed blue phase of evolu- 
tion are unknown, hence the BSG wind parameters are 
unconstrained. Even with accurate stellar parameters, 
blue loop evolution is very sensitive to assumed physics 



in the models (e.g. INeilson et alj|20lfbh . The ISM den- 
sity and the ambient magnetic field (not considered here) 
also affect shock properties, and even Betelgeuse's spa- 
tial velo city through its lo cal ISM has significant uncer- 
tainties (]Ueta et al.l 120081 ). Stronger observational con- 
straints on the ambient ISM properties would greatly 
help in constructing better models tailored specifically 
to Betelgeuse. Multiwavelength observations of the bar 
and bow shock would better constrain their dust proper- 
ties, providing more reliable mass estimates. We expect 
that new IR observations from e.g. Herschel and SOFIA 
will strongly test and constrain our model. 

5. CONCLUSIONS 

We have presented the first quantitative calculation of 
bow shock evolution for a massive runaway star undergo- 
ing the transition from blue to red on the HR diagram, 
here for a 15 M Q star evolving from BSG to RSG. In 
contrast to the previously-studied case of the transition 
from RSG to WR, here the initial wind is faster than the 
following wind and so the bow shock has time to stall and 
collapse in on itself before it is destroyed and re-formed 
by the newly expanding RSG wind. 

The BSG reverse shock collapses in on the RSG wind 
very rapidly (at the hot-bubble sound speed) , producing 
an inner shocked-shell around the expanding RSG wind, 
and eventually wraps around to completely enclose it. 
There is a short timeframe (10 kyr in this calculation) 
during which this inner shocked-shell resembles a bow 
shock and can be observed simultaneously with the much 
larger-scale bow shock left over from the BSG wind. This 
bow shock is then overtaken by the slowly expanding 
RSG wind. The transition time during which the BSG 
bow shock is still interacting with the expanding RSG 
wind can be 50 — 100 kyr, depending on the star's veloc- 
ity through the ISM, the strength of the fast wind, and 
the ambient ISM density. This is already a significant 
percentage of the time a massive star spends as a RSG, 
so it is possible that runaway RSGs may be commonly 
observed to have such features in their CSM. 

The inner shocked-shell's mass, shape, thickness and 
size are consistent with the observed properties of Betel- 
geuse's bow shock, in contrast to steady-state bow shocks 
which are typically more than an order of magnitude 
more massive. Our model also naturally provides a possi- 
ble counterpart to the infrared bar ahead of Betelgeuse's 
bow shock, the BSG wind's shocked-shell. We have also 
shown that dust emission from the inner shocked-shell 
and outer bow shock may have similar emissivities. 
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